package com.marco;

import java.util.HashMap;
import java.util.Map;

/**
 * com.marco.GisUtils.java
 * Create on 2017/3/22
 * Copyright(c) JD Inc 2017. All Rights Reserved.
 * <pre>
 * 功能说明：在此封装处理GPS 数据相关的方法
 * 1.此封装类中包含处理坐标转换的方法，使用的坐标转换算法为高斯-克吕格投影，椭球系指定坐标为WGS84坐标系
 * </pre>
 */
public class GisUtils {

    //WGS84地球椭球长半轴
    private static final Double a = 6378137.00;
    //WGS84地球椭球短半轴
    private static final Double b = 6356752.31;
    //WGS84椭球第一离心率的平方 (Math.pow(a,2)-Math.pow(b,2))/Math.pow(a,2)
    private static final Double e1 = 0.00669437999013;
    //WGS84椭球第二离心率的平方(Math.pow(a,2)-Math.pow(b,2))/Math.pow(b,2)
    private static final Double e2 = 0.006739496742227;
    //地球平均半径
    private static final Double r = 6371001.00;
    //圆周率
    private static final Double PI = 3.1415926535898;


    /**
     * 根据指定的中央子午线，进行高斯克吕格投影
     *
     * @param lon      需要正算的经度
     * @param lat      需要正算的纬度
     * @param primeLon 指定的中央子午线经度，若为空，则默认使用6度带投影
     * @return 返回对应的平面坐标系，Y值存在负值，在此不进行正值处理
     */
    public static Map<String, Double> toXY(Double lon, Double lat, Integer primeLon) {
        if (primeLon == null) {
            //采用6度带【从格林尼治子午线开始自西向东，每隔6度为1带 6N-3】投影
            primeLon = 6 * getCentralNo(lon, null) - 3;
        }
        Map<String, Double> res = new HashMap<>();
        //弧度纬度
        Double B = toRad(lat);
        //弧度经差
        Double l = toRad(lon) - primeLon * PI / 180D;
        //卯酉圈曲率半径
        Double N = getPrimeVertical(B);
        //子午线弧长
        Double lB = getMeridianLen(B);
        //常量t
        Double t = Math.tan(B);
        //常量n的平方
        Double n2 = e2 * Math.pow(Math.cos(B), 2);
        //关于纬度公式的公共部分
        Double fB = Math.cos(B) * l;
        //计算X
        Double X = lB + N * t * (Math.pow(fB, 2) / 2 + Math.pow(fB, 4) * (5 - t * t + 9 * n2 + 4 * n2 * n2) / 24
                + Math.pow(fB, 6) * (61 - 58 * t * t + t * t * t * t) / 720);
        //计算Y
        Double Y = N * (fB + Math.pow(fB, 3) * (1 - t * t + n2) / 6
                + Math.pow(fB, 5) * (5 - 18 * t * t + Math.pow(t, 4) + 14 * n2 - 58 * n2 * t * t) / 120);
        res.put("X", X);
        res.put("Y", Y);
        return res;
    }


    /**
     * 根据给定的平面坐标系，使用高斯克吕格投影进行反算到WGS84椭球坐标系
     *
     * @param x        平面坐标X
     * @param y        平面坐标Y 注意此处Y为投影带的真实值，不包含正值处理
     * @param primeLon 指定的中央子午线
     * @return WGS84椭球坐标系的经纬度
     */
    public static Map<String, Double> toLB(Double x, Double y, Integer primeLon) {
        Map<String, Double> res = new HashMap<>();
        //获取底点纬度，以大地测量角度计算极限误差作为迭代限差
        Double Bf = getBottomLat(x, toRad(0.000001));
        //常数t
        Double t = Math.tan(Bf);
        //常数n的平方
        Double n2 = e2 * Math.pow(Math.cos(Bf), 2);
        //常数V的平方
        Double v2 = 1 + e2 * Math.pow(Math.cos(Bf), 2);
        //卯酉圈半径
        Double N = getPrimeVertical(Bf);

        //计算纬度
        Double B = Bf - v2 * t * (Math.pow(y / N, 2) - Math.pow(y / N, 4) * (5 + 3 * t * t + n2 - 9 * n2 * t * t) / 12
                + Math.pow(y / N, 6) * (61 + 90 * t * t + 54 * t * t * t * t) / 360) / 2;

        //计算经度
        Double l = (y / N - Math.pow(y / N, 3) * (1 + 2 * t * t + n2) / 6 +
                Math.pow(y / N, 5) * (5 + 28 * t * t + 24 * Math.pow(t, 4) + 6 * n2 + 8 * t * t * n2) / 120)
                / Math.cos(Bf);
        B = toAngle(B);
        l = toAngle(l) + primeLon;
        res.put("lat", B);
        res.put("lon", l);
        return res;
    }


    /**
     * 获取中央子午线的投影带带号
     *
     * @param lon  经度
     * @param wide 指定带宽，若为空，默认使用6度带
     * @return 高斯投影带带号
     */
    private static Integer getCentralNo(Double lon, Integer wide) {
        if (wide == null) {
            wide = 6;
        }
        Integer index = 0;
        switch (wide) {
            case 3:
                index = (int) ((lon - 1.5) / wide);
                break;
            case 6:
                index = (int) (lon / wide);
                break;
        }
        return index + 1;
    }


    /**
     * 利用牛顿迭代算法获取指定X坐标的底点纬度，
     *
     * @param x   平面坐标X
     * @param eps 指定限差，单位为弧度
     * @return 返回底点纬度
     */
    private static Double getBottomLat(Double x, Double eps) {
        Double m0 = r * (1 - e1);
        Double m2 = 3 * e1 * m0 / 2D;
        Double m4 = 5 * e1 * m2 / 4D;
        Double m6 = 7 * e1 * m4 / 6D;
        Double m8 = 9 * e1 * m6 / 8D;

        Double a0 = m0 + m2 / 2 + 3 * m4 / 8 + 5 * m6 / 16 + 35 * m8 / 128;
        Double a2 = m2 / 2 + m4 / 2 + 15 * m6 / 32 + 7 * m8 / 16;
        Double a4 = m4 / 8 + 3 * m6 / 16 + 7 * m8 / 32;
        Double a6 = m6 / 32 + m8 / 16;
        Double a8 = m8 / 128;

        Double preB = 0D;
        Double nextB = x / a0;
        while (Math.abs(nextB - preB) > eps) {
            preB = nextB;
            nextB = (x + a2 * Math.sin(2 * preB) / 2 - a4 * Math.sin(4 * preB) / 4 + a6 * Math.sin(6 * preB) / 6
                    - a8 * Math.sin(8 * preB) / 8) / a0;
        }
        return nextB;
    }


    /**
     * 角度转弧度
     *
     * @param angle 角度
     * @return 弧度
     */
    private static Double toRad(Double angle) {
        return angle * PI / 180D;
    }

    /**
     * 弧度转角度
     *
     * @param rad 弧度
     * @return 角度十分制
     */
    private static Double toAngle(Double rad) {
        return rad * 180 / PI;
    }


    /**
     * 根据制定的纬度（弧度制）计算卯酉圈曲率半径
     *
     * @param latRad 弧度纬度
     * @return 卯酉圈半径
     */
    private static Double getPrimeVertical(Double latRad) {
        Double n0 = r;
        Double n2 = e1 * n0 / 2D;
        Double n4 = e1 * n2 * 3 / 4D;
        Double n6 = e1 * n4 * 5 / 6D;
        Double n8 = e1 * n6 * 7 / 8D;
        Double sinB = Math.sin(latRad);
        //推导公式结果
        return n0 + n2 * Math.pow(sinB, 2) + n4 * Math.pow(sinB, 4) + n6 * Math.pow(sinB, 6) + n8 * Math.pow
                (sinB, 8);
    }

    /**
     * 根据给定的纬度计算子午线弧长
     *
     * @param latRad 弧度制的纬度
     * @return 子午线弧长
     */
    private static Double getMeridianLen(Double latRad) {
        Double m0 = r * (1 - e1);
        Double m2 = 3 * e1 * m0 / 2D;
        Double m4 = 5 * e1 * m2 / 4D;
        Double m6 = 7 * e1 * m4 / 6D;
        Double m8 = 9 * e1 * m6 / 8D;

        Double a0 = m0 + m2 / 2 + 3 * m4 / 8 + 5 * m6 / 16 + 35 * m8 / 128;
        Double a2 = m2 / 2 + m4 / 2 + 15 * m6 / 32 + 7 * m8 / 16;
        Double a4 = m4 / 8 + 3 * m6 / 16 + 7 * m8 / 32;
        Double a6 = m6 / 32 + m8 / 16;
        Double a8 = m8 / 128;
        //推导公式结果
        return a0 * latRad - a2 * Math.sin(2 * latRad) / 2 + a4 * Math.sin(4 * latRad) / 4
                - a6 * Math.sin(6 * latRad) / 6 + a8 * Math.sin(8 * latRad) / 8;
    }


}
